Wall-Enhanced Convection in Vibrofluidized Granular Systems 
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An event-driven molecular dynamics simulation of inelastic hard spheres contained in a cylin- 
der and subject to strong vibration reproduces accurately experimental results jl| for a system of 
vibrofluidized glass beads. In particular, we are able to obtain the velocity field and the density 
and temperature profiles observed experimentally. In addition, we show that the appearance of 
convection rolls is strongly influenced by the value of the sidewall-particle restitution coefficient. 
Suggestions for observing more complex convection patterns are proposed. 



PACS numbers: 45.70.Mg, 47.20.Bp, 47.27.Te, 81.05.Rm 
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If sufficient power is supplied to a granular system con- 
tained in a cylindrical container by a vertically vibrating 
wall, fluidization occurs. Once a steady state is obtained, 
the energy input to the system by the wall is, on average, 
balanced by the energy dissipated in inelastic grain-grain 
and grain- wall collisions. This energy flux is responsi- 
ble for various effects including non-linear temperature 
profiles, heap formation^ and convection ||f|, ^. Re- 
cent experiments using positron emission particle track- 
ing (PEPT) |]l] have observed buoyancy-driven convection 
in a highly fluidized granular system for a range of grain 
numbers and shaker amplitudes. While convection effects 
have been observed in a two-dimensional model system 
[Q, 1^ , no simulation studies of a three dimensional system 
corresponding to the PEPT experiments have yet been 
reported. 

We study the behavior of a model system in which the 
grains are modeled as inelastic hard spheres. Most of the 
previous studies of this model focused on the homoge- 
neous cooling state, clustering, and kinetic and hydrody- 
namic theory for particles with small inelasticity 0, ||. 
Unlike many other many-body systems, the model stud- 
ied here closely matches the experimental system. In 
particular, the number of particles is identical and we 
use values of particle-particle and particle-wall restitu- 
tion coefficients that were determined by independent 
experimental measurements. We show that the model 
reproduces accurately the experimentally observed den- 
sity and granular temperature profiles||^. Moreover, for 
a range of conditions, a velocity field exists for which 
there is a net circulation. Interestingly, the direction and 
intensity of the circulation are strong functions of the 
particle-wall restitution coefficient. In addition, we sug- 
gest a modification to the experiment that may generate 
new convection patterns. 



The system consists of a number N of hard spheres 
contained in a cylinder of radius R. The spheres col- 
lide inelastically with each other and with the side walls 
with coefficients of restitution c and , respectively. Be- 
tween collisions, the spheres are subject to a downward 
constant acceleration due to the (vertical) gravitational 
field. Energy is injected into the system by the bottom 
wall which vibrates with a symmetric saw-tooth profile 
characterized by an amplitude A and a period T. We 
do not expect the behavior of the system to be strongly 
dependent on the form of the profile jiot. 

The equations describing the dynamics, which follow 
from conservation of momentum, for the particle side- 
wall, particle bottom-wall and particle particle collisions 



are: 
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-[(vj - Vi).n]n 



(1) 



(2) 



(3) 



respectively where v^^^ is the velocity of particle i in the 
X, y plane, w„ is the velocity of the bottom wall, Vi^z is 
the z vertical component of the velocity of particle i, f 
is the unit position vector of particle i and n is the unit 
center-to-center vector between the colliding pair i and j. 
Note that we have taken the tangential restitution coeffi- 
cient for sphere-sphere, and sphere-wall, collisions equal 
to one. The normal restitution coefficient for sphere- 
bottom wall collisions was also taken as unity. 

One problem with the event-driven simulation is that 
a sphere may collide more and more frequently with the 
side wall as the radial component of its momentum is 
dissipated. To avoid the inelastic collapse [Q, which has 
been observed in simulations of bulk granular systems. 
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FIG. 1: Mean velocity field for a vibrofiuidized system of 
inelastic hard spheres in the (r, z) plane. The wall-particle 
restitution coefficient is c„ = 0.68. Other parameters are 
given in the text 

we simply inject a small amount of energy in the radial 
direction when the radial velocity falls below a cut-off 
value. This condition arises infrequently and its handling 
in this way has no discernible effect on the results. 

To reproduce the experimental conditions of Wildman 
et al j|] we used parameter values TV — 1050, c = 0.91, 
Cw = 0.68. Taking the unit of length as the sphere di- 
ameter, (T, the cylinder radius is i?/cr = 14. Since the 
base-wall particle coefficient is not given, we assume a 
value of one. The shaker amplitude is thus left as the 
only adjustable parameter for the comparison of the sim- 
ulation with the experiments. We performed different 
runs with a dimensionless shaking amplitude varying be- 
tween 0.05 and 0.6. To generate the initial configuration, 
spheres were inserted sequentially and randomly in the 
cylinder without overlap. A preliminary simulation was 
then performed, typically for 5000 collisions per parti- 
cle, in order to allow the system to reach the stationary 
non-equilibrium state. 

Figure |l| displays the mean velocity field for the sys- 
tem of inelastic spheres in the (r, z) plane (r denotes 
the radial distance from the axis of the cylinder) when 
the reduced shaker amplitude is A/a = 0.32. The field 
is averaged over the azimuthal angle and approximately 
2 X 10^ collisions). A torroidal convection roll is clearly 
present in which the particles flow, on average, up from 
the center and down the side wall. Note that, unlike the 
related Rayleigh-Benard phenomenon, the roll is asym- 
metric. The velocity field is in near quantitative agree- 
ment with that observed experimentally by Wildman et 
al. 01 . In particular, the center of the vortex roll is at 
(r/(T ~ 10, z/cr ~ 6) in reduced units which corresponds 



FIG. 2: Same as Figure |l| except that the wall-particle resti- 
tution coefficient is Cw = 1.0. Note that the velocities are 
much smaller than in Figure ^. 

to (r ~ 50mm, z 2± 30mm) in experiments [Q. The shape 
of the roll, as well the location of its center, do not change 
appreciably as the shaker amplitude is varied between 
0.28 < A/a < 0.37. 

To highlight the role played by the wall, we performed 
additional simulations in which the wall-sphere collision 
is taken as elastic (a "perfect' wall for which the wall- 
particle restitution coefficient is equal to 1) with the other 
parameters unchanged. The mean velocity field of this 
simulation is quite different from the inelastic case. Al- 
though some convection is apparent in Figure |^, (parti- 
cles on average flow up the wall and down to the center). 
It is weaker by about an order of magnitude in intensity 
(compare the velocity units at the top of Figs. |l] and. H) 
and the roll direction is opposite. 

To quantify the bulk convection, we have calculated 
the total velocity correlation, Ctot — {^/^c)'^i,jC{i,j)^ 
where (i,j) label the cells in the hydrodynamic region, 
Nc their total number and 

C(*,j) = (l/8)I],,,,vv(z,j).v(j',/) (4) 

where the sum is over the 8 neighboring cells of 

(i,j). See Figure^. Bulk convection is present when Ctot 
is different from zero. 

For Cw = 0.68, as the shaker amplitude is increased 
the total velocity correlation Ctot becomes different from 
zero when A/ a > 0.05, increases until a/a ~ 0.2 and 
stays roughly constant for larger amplitudes. A similar 
behavior was observed in the experiments of Wildman et 
al. jl]. For Cm = 1, Ctot is almost equal to zero for all 
amplitudes (Figure ^ . A close examination of the corre- 
sponding flow velocity flelds show that the amplitude of 
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FIG. 3: Order parameter Ctot versus the shaker amphtude 
A/a oi the base waU. The squares and circles correspond to 
simulations where the wall-particle restitution coefficient is 
Cu, = 1.0 and 0.68, respectively . 



the circulation of particles up the wall and down to the 
center exists for A/ a > 0.05, a feature which is not well 
handled by the calculation of the order parameter. 

Figure ^ shows the granular temperature as a function 
of the altitude, z, for different values the radial posi- 
tion. These results agree well with the experiments (Fig.l 
of Ref^) except for an underestimated maximum close 
to the basewall. These differences could result from the 
base-particle collisions which are assumed elastic in sim- 
ulations. A minimum in the granular temperature profile 
is observed in the current simulations (and also in exper- 
iments) in the neighborhood of the center of the cylinder, 
but becomes weaker close to the wall cylinder. Fig |] dis- 
plays the same profiles for the system with c^, = 1: The 
temperature minimum in the vertical direction is more 
pronounced, but the temperature varies little in the ra- 
dial direction. The existence of a minimum in the granu- 
lar temperature along the z-direction is a phenomenon 
which is also present in the absence of convection 
and results from the fact that at high altitudes the den- 
sity and collision frequency are small, leading to a small 
number of hot particles. 

Figure |^ displays a packing fraction plot which agrees 
qualitatively with the experimental results, except in the 
bottom region of the side wall where the density is higher 
than in experiments. Nevertheless, there is a maximum 
in the density profile in the vertical direction, whatever 
the radial distance. Moreover, there is also a maximum 
along the radial direction which is a consequence both 
of the compressible nature of the inelastic hard sphere 




FIG. 4: Granular temperature for a vibrofluidized sys- 
tem of inelastic hard spheres as a function of the al- 
titude z for different values of the radial position. 
From bottom to top, the curves correspond to r/a = 
1.3, 2.6, 3.9, 5.2, 6.5, 7.8, 9.1, 10.4, 11.7, 13. The particle mass 
is 1.875 X 10"''%. 

system and of the small wall-particle restitution coeffi- 
cient. When = 1, the radial density gradient almost 
vanishes. 

Although a complete theoretical treatment of the con- 
vection is beyond the scope of this Letter, we present 
some simple arguments that may help to rationalize 
the observed behavior. The Rayleigh-Benard convection 
which is the paradigm for pattern formation in classi- 
cal fluid mechanics, has been extensively studied over 
the past three decades p^, p^ . A temperature gradient 
applied to a fluid leads to a competition between buoy- 
ancy and disspation, which leads to the occurence of rolls 
above a threshold characterized by the non dimensional 
Rayleigh number Ra , which characterizes the ratio of the 
two forces. For larger temperature gradients, on observes 
a cascade of bifurcations ending in a chaotic regime. 

For classical fluids, hydrodynamic equations reproduce 
the observed patterns remarkably well. The range of va- 
lidity of continuum equations applied to granular materi- 
als, which contain a sink term in the energy equation to 
account for the dissipation of energy during collisions, is 
more restricted. Hydrodynamic equations can provide a 
good description of inelastic spheres ||l^ , provided that 
the particles packing fraction remains 77 > 0.04 (below 
which the mean free path becomes comparable to the 
diameter of the cylinder). Similarly, if the packing frac- 
tion is too large (ry > 0.3) the dissipation becomes too 
strong and the description breaks down. Using such an 
approach, Brey et al.[0 obtained temperature profiles 
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FIG. 5: Same as Figure ^ except Cw = 1.0. From 
top to bottom, the curves correspond to r/a = 
13, 11.7, 10.4, 9.1, 7.86.5, 5.2, 3.9, 2.6, 1.3. 
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FIG. 6: Packing fraction (77) plot for the vibrofiuidized sys- 
tem of inelastic hard spheres, with Cw — 0.68. The greyscale 
step corresponds to an rj increase of 0.02 and white regions 
correspond to r; < 0.02. 

for an open granular system, without convection and pre- 
dicted the presence of a minimum as a function of 2, 
coinciding with the limit of validity of the hydrodynamic 
equations. 

When Cw = 0.68, the hydrodynamic description breaks 
down close to the wall (high density and temperature gra- 
dient, strong dissipation). The onset of the bulk convec- 
tive regime with a "real" wall is lowered because of the 
existence of the radial temperature gradient in addition 



to the vertical one. This can explain both the direction 
of circulation and the lowering of the threshold of convec- 
tion compared to the case when the wall is perfect and no 
bulk convection is present. The circulation of particles is 
restricted close to the sidewall where strong dissipation is 
not well described by hydrodynamic equations. A kinetic 
approach is required close to the cylinder. 

Our simulations clearly demonstrate the relevance of 
the inelastic hard sphere model for describing the exper- 
imental study of Wildman et al[Q and underline the key 
influence of the cylinder wall on the convection. Experi- 
mentally, unlike in the simulations, the wall-particle resti- 
tution coefficient is not an easily tunable parameter. We 
therefore propose a modification to the experiment that 
may lead to the observation of new phenomena. Specif- 
ically, we suggest adding an inner (solid) cylinder of a 
given radius, coaxial with the original cylinder. In the 
neighborhood of this inner cylinder, we expect particles 
to have a downward motion, and two, or more, toroidal 
rolls with alternate directions of circulation, as in a fluid 
system, might be produced. 
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